# Title     : TODO
# Objective : TODO
# Created by: Administrator
# Created on: 2019/6/18

library(gplots)
library(Hmisc)
library(optparse)
library(dplyr)
library(ropls)
library(magrittr)
library(tibble)
library(stringr)
library(ComplexHeatmap)
library(tibble)
library(tidyr)
library(tools)

option_list <- list(
make_option("--i", default = "AllMet.csv", type = "character", help = "metabolite data file"),
make_option("--g", default = "true_SampleInfo.csv", type = "character", help = "sample group file"),
make_option("--et", default = "", type = "character", help = "extra data dir"),
make_option("--confounder", default = "", type = "character", help = "confounder data file")
)
opt <- parse_args(OptionParser(option_list = option_list))

createWhenNoExist <- function(f){
    ! dir.exists(f) && dir.create(f)
}

zd_pcor <- function (left = left , right =right  , confounding = confounding , unique_id_colname = "ID" , mannual_n_for_fdr = F, filtered_pairs_saved_in = "global_filtered_pairs"  ,
method = "spearman"
){
    if (exists(filtered_pairs_saved_in) == F) {
        global_filtered_pairs <<- data.frame()
        cat ("初始化了相应的变量,将剔除掉的变量对子保存在了 global_filtered_pairs 变量中. 这个名字是固定写死的,不要更改.")
    }

    part_1 <- left[, unique_id_colname] %>% as.character()
    part_2 <- right[, unique_id_colname] %>% as.character()
    part_3 <- confounding[, unique_id_colname] %>% as.character ()

    pooled_ID <- c(part_1, part_2, part_3) %>% unique()
    neworder <- pooled_ID %>% order()
    pooled_ID <- pooled_ID[neworder] %>% data.frame(. , stringsAsFactors = F)
    colnames(pooled_ID) <- unique_id_colname

    left <- merge (pooled_ID , left , by = c(unique_id_colname) , all.x = T)
    right <- merge (pooled_ID , right , by = c(unique_id_colname) , all.x = T)
    confounding <- merge (pooled_ID , confounding , by = c(unique_id_colname) , all.x = T)

    left <- left[, - 1]
    right <- right[, - 1]
    confounding <- confounding[, - 1]

    for (i in 1 : ncol (left)) {
        if (i == 1) { pooled_r_and_p <- data.frame()
            pooled_error <- data.frame ()
        }
        for (j in 1 : ncol (right)) {
            pair_1 <- colnames (left)[i]
            pair_2 <- colnames (right)[j]

            temp_frame <- cbind (left[, i], right[, j], confounding)
            temp_frame <- na.omit(temp_frame)
            temp_frame <- apply (temp_frame , 2, function (each_col){
                result <- each_col %>% as.character () %>% as.numeric()
            })

            IS_error <<- F
            if (T) {
                tryCatch(each_line <- ppcor::pcor.test(x = temp_frame[, 1], y = temp_frame[, 2], z = temp_frame[, 3 : ncol(temp_frame)], method = c(method)),
                error = function(e) {
                    IS_error <<- T
                    cat (i , ' ', pair_1 , j, ' ', pair_2 , "的偏相关系数无法被计算\n")

                    each_filtered_pair <- data.frame (pair_1 = pair_1, pair_2 = pair_2)
                    global_filtered_pairs <<- rbind (global_filtered_pairs, each_filtered_pair)
                })
            }

            if (IS_error == F) {
                each_line <- data.frame (pair_1 = pair_1 , pair_2 = pair_2 , each_line , stringsAsFactors = F)
                colnames(each_line) <- c("pair_1", "pair_2", "r", "p", "statistic", "n", "gp", "Method")
                pooled_r_and_p <- rbind (pooled_r_and_p , each_line)
            }
        }
    }

    IS_filtered <- is.na (pooled_r_and_p[, "p"]) == T
    temp_filtered_data <- pooled_r_and_p[IS_filtered  , c("pair_1", "pair_2")]
    global_filtered_pairs <<- rbind (global_filtered_pairs , temp_filtered_data)
    pooled_r_and_p <- pooled_r_and_p[! IS_filtered,]

    if (mannual_n_for_fdr > 0) {n <- mannual_n_for_fdr} else {
        n <- nrow(pooled_r_and_p)}
    pooled_r_and_p$fdr = p.adjust (pooled_r_and_p$p, method = "fdr" , n = n)
    pooled_r_and_p$BH = p.adjust (pooled_r_and_p$p, method = "holm", n = n)

    return (pooled_r_and_p)
}

sampleIds <- read.csv(opt$g, header = T, stringsAsFactors = F) %>%
    select(c("SampleID", "ClassNote")) %>%
    .$SampleID

diffData <- read.csv("../../07/Sig_Met.csv", header = T, stringsAsFactors = F, row.names = 1) %>%
    rownames_to_column("Metabolite") %>%
    filter(IS_Final_Pooled_Sig == 1)

diffNames <- diffData %>%
.$Name

data <- read.csv(opt$i, header = T, stringsAsFactors = F, check.names = F) %>%
    select(- c("HMDB", "KEGG", "Class")) %>%
    select(c("Metabolite", sampleIds)) %>%
    filter(Name %in% diffNames) %>%
    column_to_rownames("Metabolite") %>%
    as.data.frame() %>%
    t()

files <- list.files(opt$et, full.names = T)
# for (file in files) {
for (i in 1:1) {
    file=files[i]
    parent <- basename(file) %>%
    file_path_sans_ext()
    createWhenNoExist(parent)
    extraData <- read.csv(file, row.names = 1, check.names = F) %>%
        rownames_to_column("Metabolite") %>%
        select("Metabolite", sampleIds) %>%
        column_to_rownames("Metabolite") %>%
        t() %>%
        as.data.frame() %>%
        rownames_to_column("ID")
    print(head(extraData))
    confounderData <- read.csv(opt$confounder, check.names = F,row.names=1)# %>%
        # rownames_to_column("Metabolite") %>%
        # select("Metabolite", sampleIds) %>%
        # column_to_rownames("Metabolite") %>%
        # t()
    print(head(confounderData))
    listRs <- zd_pcor(data, extraData, confounderData)

    allData <- data.frame(node1 = listRs$node1, node2 = listRs$node2, cor = listRs$cor, p = listRs$p,
    stringsAsFactors = F) %>%
        mutate_at(vars(c("cor")), function(x){
            ifelse(is.na(x), 0, x)
        }) %>%
        mutate_at(vars(c("p")), function(x){
            ifelse(is.na(x), 1, x)
        }) %>%
        mutate(fdr = p.adjust(p, method = "fdr"))

    write.csv(allData, paste0(parent, "/for_Cytoscape.csv"), quote = T, row.names = F)

    corData <- allData %>%
        select(c("node1", "Variable", "cor")) %>%
        spread(node1, "cor")
    write.csv(corData, paste0(parent, "/r_values_all.csv"), quote = T, row.names = F)

    pData <- allData %>%
        select(c("node1", "Variable", "p")) %>%
        spread(node1, "p")
    write.csv(pData, paste0(parent, "/p_values_all.csv"), quote = T, row.names = F)

    fdrData <- allData %>%
        select(c("node1", "Variable", "fdr")) %>%
        spread(node1, "fdr")
    write.csv(fdrData, paste0(parent, "/fdr_values_all.csv"), quote = T, row.names 